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Abstract 

Background: The CRISPR-Cas systems of adaptive antivirus immunity are present in most arcliaea and many 
bacteria, and provide resistance to specific viruses or plasmids by inserting fragments of foreign DNA into tlie liost 
genome and tlien utilizing transcripts of tinese spacers to inactivate tlie cognate foreign genome. Tlie recent 
development of powerful genome engineering tools on the basis of CRISPR-Cas has sharply increased the interest in 
the diversity and evolution of these systems. Comparative genomic data indicate that during evolution of prokaryotes 
CRISPR-Cas loci are lost and acquired via horizontal gene transfer at high rates. Mathematical modeling and initial 
experimental studies of CRISPR-carrying microbes and viruses reveal complex coevolutionary dynamics. 

Results: We performed a bifurcation analysis of models of coevolution of viruses and microbial host that possess 
CRISPR-Cas hereditary adaptive immunity systems. The analyzed Malthusian and logistic models display complex, 
and in particular, quasi-chaotic oscillation regimes that have not been previously observed experimentally or in 
agent-based models of the CRISPR-mediated immunity. The key factors for the appearance of the quasi-chaotic 
oscillations are the non-linear dependence of the host immunity on the virus load and the partitioning of the 
hosts into the immune and susceptible populations, so that the system consists of three components. 

Conclusions: Bifurcation analysis of CRISPR-host coevolution model predicts complex regimes including 
quasi-chaotic oscillations. The quasi-chaotic regimes of virus-host coevolution are likely to be biologically 
relevant given the evolutionary instability of the CRISPR-Cas loci revealed by comparative genomics. The results 
of this analysis might have implications beyond the CRISPR-Cas systems, i.e. could describe the behavior of any 
adaptive immunity system with a heritable component, be it genetic or epigenetic. These predictions are 
experimentally testable. 

Reviewers' reports: This manuscript was reviewed by Sandor Pongor, Sergei IVlaslov and Marek Kimmel. For the 
complete reports, go to the Reviewers' Reports section. 



Background 

The arms races between microbes and viruses preying on 
them often display rich, complex population dynamics [1]. 
In principle, the dynamics of virus-microbe interactions is 
analogous to the classical predator-prey models [2-5] 
but both microbes and viruses evolve much faster than 
animals such that virus-host interactions change on a 
scale that may be amenable to direct laboratory study. 
One of the adaptation mechanisms employed by hosts 
to curb viruses is the CRISPR-Cas (Clustered 7?egularly 
interspaced Short Palindromic 7?epeats-CRISPR associated 
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proteins), a recentiy discovered adaptive immunity system 
that is present in the great majority of Archaea and 
many bacteria [6-12]. Microbes create heritable memory 
of viruses that attack them by inserting virus-derived 
spacers into CRISPR repeat cassettes, thus following the 
Lamarckian modality of evolution that dramatically accel- 
erates adaptation [13]. The rapid adaptation through the 
activity of CRISPR-Cas is possible because this system en- 
genders heritable genetic changes that are directly benefi- 
cial for the archaeon or bacterium in the face of a specific 
environmental challenge (a virus), in contrast to the ran- 
dom, undirected mutations in the Darwinian evolutionary 
framework [14]. The CRISPR-Cas systems are increasingly 
used as powerful, versatile tools for genomic engineering 
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tools which sharply increases the interest in their diversity 
and evolution [15-20]. 

General considerations suggest that the population 
dynamics of virus-host coevolution should be dominated 
by periodic selective sweeps alternating between the host 
(when it "discovers" a resistance mutation or acquires 
immunity against the dominant virus lineage) and the 
virus (when an immunity escape mutation occurs), similar 
to the case of rapidly evolving human viruses [21]. Indeed, 
such behavior has been observed in simple Lotka-Volterra 
type models of phage-bacteria coevolution [22,23] as well 
as in direct evolutionary experiments [23] . However, direct 
and indirect population studies reveal much more com- 
plex behaviors of the actual populations that usually does 
not involve strain dominance and instead displays long- 
term persistence of multiple lineages of both the microbial 
hosts and the viruses [24-27]. 

Several detailed agent-based models of coevolution 
between viruses and CRISPR-Cas-carrying hosts have 
been developed and analyzed [26,28-35]. The agent-based 
models incorporate the salient features of the CRISPR-Cas 
system such as the existence of the CRISPR cassette with 
virus-derived spacers, immunity conferred to a host by 
spacers that match the attacking virus and acquisition of 
new spacers as a result of failed virus infections. These 
models allow one to reproduce many aspects of the 
observed behavior of coevolving virus-host systems and 
predict conditions required for the evolutionary mainten- 
ance of the CRISPR-Cas immunity, such as a threshold of 
viral diversity [28]. 

Agent-based models provide for the exploration of 
interactions of arbitrary complexity and naturally incorp- 
orate the desired level of granularity (e.g. individual-based 
or lineage-based models) and the stochasticity of the 
processes involved. However, such models typically 
possess a high-dimensional parameter space that cannot 
be explored in full, so that not all potential regimes, some 
of which could be biologically relevant, are captured. 
In contrast, mathematical models based on systems of 
differential equations are limited in complexity and are 
inherently less realistic (at the very least because they 
approximate reality with infinitely small deterministic 
changes) but when analytically tractable, permit a full 
and rigorous analysis of all possible behaviors. 

Here we describe Lotka-Volterra type models of 
interaction between a host with a heritable adaptive 
immunity system, such as CRISPR-Cas, and a virus 
that escapes the immunity via implicit accumulation of 
mutations which is implemented as gradual immunity 
decay. We construct "minimal" analytical models which 
capture qualitatively the basic regimes of the CRISPR-Cas 
system behaviors previously found experimentally and 
through the agent-based modeling. We explore the full 
spectrum of possible behaviors of this virus-host system. 



compare the results with those of a more detailed agent- 
based model [32], and describe a previously unnoticed 
regime of quasi-chaotic oscillations. 

Results 

Three-component CRISPR population dynamics: 
Malthusian and logistic versions 

In order to construct "minimal" analytical models that 
would qualitatively capture the basic regimes of the 
CRISPR-Cas system behaviors we analyzed, as a prelim- 
inary step, a two-component Volterra-type model which 
describes the dynamics of virus-host system and con- 
sider immunity static within the timescale of the model. 
Our analysis shows that two-component models display 
simple dynamical regimes (see Methods for details). 
However, the CRISPR-Cas system of adaptive immunity, 
which this work seeks to model, cannot be expected to 
follow these simple regimes given its dynamic evolution 
that involves rapid acquisition of immunity to a particu- 
lar virus or plasmid along with loss and gain of entire 
CRISPR-Cas loci [36]. Thus, more realistic models 
should take into account that adaptive immunity sys- 
tems are characterized by the existence of immune 
memory that enhances the response in individuals en- 
countering a familiar challenge [37,38]. Originally non- 
immune ("native") individuals can acquire the adaptive 
response capacity that persists at timescales relevant for 
our model. Thus, we introduce two categories of hosts, 
immune and non-immune, that interact differently 
with the virus and exchange individuals via immunity 
acquisition and decay. 

Let x{t) be the density of immune hosts with the im- 
munity p with respect to viruses, 0 <p <1; y{t) be the 
density of sensitive hosts with immunity s, 0 <s <p<l; 
z{t) be the density of viruses. We consider the following 
3-component model: 

dx 

— = x{\-l-a{x -\- y))-bxz{\-p) + esyz = P{x,y,z), 
at 

dy 

— = 7 + lx-ay{x + y)-byz{l-s)-esyz = Q{x,y,z), 
dz 

— = z(-d + bM{x{l-p) + y(l-s))-b{xp + ys)) =R(x,y,z). 
dt 

(1) 

Here / is the immunity decay rate (x — > j flow), e is the 
immunity acquisition rate; d is the death rate of viruses, 
M is the virus reproduction rate; b is the encounter rate 
coefficient; the growth rates of immune and sensitive 
hosts are equal to 1. 

Model (1) for a = Q describes a situation when immune 
and sensitive hosts in the absence of viruses grow according 
to the Malthusian model. If a > 0, the model (1) describes a 
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more realistic situation when both classes of hosts grow 
according to the logistics model. 
Coordinates of equilibria of (1) solve the system: 

P{x,y,z) = x{l-l-a{x + y))-bxz{l-p) + esyz = 0, 
Q{x,y,z) =y + lx-ay{x + y)-byz{l-s)-esyz = 0, 
R(x,y,z) =z{-d + bM(x(l-p) + y(l-s))-b(xp + ys)) = 0. 

(2) 

Model (1) has equilibrium 0(0, 0, 0) in all cases. The 
logistic model with a> 0 has one more equilibrium A 
(0,l/fl,0), which corresponds to existence of only non- 
immune hosts. 

In Methods the following assertions are proven. 

Statement 1. 



Let us consider model (1) with the immunity p = p{z) 
defined by (3). We do not attempt a complete analysis of 
this model but rather seek to identify stable modes and 
most interesting dynamical behaviors. 

We start with the Malthusian version when a = 0. It 
has non-trivial equilibrium Bg(Xe, ye, Ze) such that x = Xe, 
y = ye are expressed via z = z^- 

-d{\-b{\-s)z) d{l-b{\-p)z) 
' b{p-s){l+M-bz)' b(p-s){l+M-bz)' ^ ' 

where z = Zg solves the equation 

(1-0 + b{-2 + 1+P + 5-es(l-/))z + b^{l-p) (5) 
X (1-s + es)z^ = 0. 



(1) Equilibrium 0{x = y = z = Q) is unstable for all 
parameter values; 

(2) If a>0 then equilibrium A(0,l/a,0) is stable for 
M < and unstable for M > fi±^ . 

If the immunity p in model (1) is a constant, then, 
additionally, the model may have either non-trivial stable 
equilibrium, or may demonstrate periodic oscillations. 
Detail description of possible behaviors of the model 
with constant p is given in Methods. 

More realistically, the immunity p is not a constant but 
depends on the density of the virus, p = p{z) (0<p< 1). 
Below we consider the latter case, in agreement with 
empirical observations and computer simulations. Specif- 
ically, it has been shown that CRISPR-Cas systems are 
(nearly) ubiquitous in archaeal and bacterial hyperthermo- 
philes are present in less than half of the available mesophile 
genomes [28,32,39,40]. Analysis of agent-based models of 
virus-host coevolution suggest that this distinction stems 
from the fact that hyperthermophUes face lower virus 
loads and diversity than mesophiles providing for higher 
efficacy of CRISPR-Cas [28]. 

Our aim is to find all stable modes of the model at 
different values of the model parameters and to describe 
the transitions from one mode to another when parameters 
vary; by other words, we want to construct the bifurcation 
diagram of the model. It is natural to suppose that p = p(z) 
monotonically decreases and tends to the immunity 5 of 
sensitive hosts at large z. From now on we consider: 

piz) = {l-s)e-'" + s (3) 

where k, s are constants, 0 < s < 1, A: > 0. Under equation 
(3), immunity is a monotonically declining function of 
the virus amount that tends to a constant, maximum 
p (maximally efficient adaptive immunity) when z tends 
to zero, and tends to 5 (no adaptive immunity, innate 
immunity only) when z tends to infinity. 



In particular, we show that for a wide domain of the 
parameter values, the model demonstrates non-periodic 
oscillation of all three variables. The following assertions 
are valid (see Methods). 

Statement 2. For a wide range of (fixed) parameters I, e, 
s,b,0 <k<l, system (1), (3) with a = 0 has only one posi- 
tive and unstable equilibrium Be{xe,ye, Ze) under condi- 
tion p{Ze) <yv^; the coordinates (x^jg, zj of this 
equilibrium satisfy (4), (5). 

Trajectories of the model show quasi-chaotic behavior 
for a broad range of the model parameters. Consider in 
detail the behavior of the model solutions depending on 
the values of the parameters /, e, and M. For / > e, typical 
trajectories starting close to the equilibrium point are 
shown at Figure 1. Initially, all variables show almost 
periodic oscillations with increasing amplitudes. Then, 
as the amplitudes become large enough, the behavior of 
the trajectories changes sharply and the oscillations 
become (quasi)-chaotic; if the initial point is far from 
the equilibrium, then the (quasi) -chaotic oscillations are 
observed from the very beginning. When I <e the behav- 
ior of the trajectories is similar. The difference is that the 
fraction of immune hosts, x, in case / < e is greater than it 
is in case l> e. Again, if the initial point is taken far from 
the equilibrium, then the (quasi) -chaotic oscillations are 
observed from the very beginning, similar to Figure 1. 
These types of behavior are observed in a wide area of 
values of the parameter M, 1 <M < 1000. Notice, that 
when M increases, the maximum values of x,y decrease 
whereas z does not depend on M. This effect does not 
seem to have a plausible biological interpretation (viruses 
cannot exist if the hosts go extinct), indicative of apparent 
limitations of the model. 

Let us consider a more realistic version of 3-component 
system (1) with logistic growth of hosts and the immunity 
p{z) given by (3). Again, we do not attempt a complete 
analysis of this model but rather seek to identify stable 
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Figure 1 Solutions {x(t),y(t),z(f)} and phase portrait of the Malthusian model (1), a = 0 with parameter values /= 0.9, e = 0.1. Other 
parameters are M= 100,d = l,b = 0.01 ,/; = 0.5, 5 = 0.1. Initial values x(0) = 0.07, y(0) =0.15, z(0) =22.25 are chosen to be close to the equilibrium. 



modes and most interesting dynamical behaviors and 
to compare them to those of the Malthusian version of 
the model. In particular, we show that the model displays 
non-trivial stable equilibria, stable oscillations, and quasi- 
chaotic oscillations of all three variables. More precisely, 
there exists a critical value of the virus birth rate Af^ at 
fixed values of all other parameters such that the system 
tends to a stable equilibrium when M < Af^, but if M > 
'i then small periodic oscillations appear due to the 
Hopf bifurcation and then, if M > > A1^\ the system tran- 
sits to a regime of (quasi)-chaotic oscillations. 

Coordinates of any non-trivial equilibrium B{xi,,ye,Ze) 
should satisfy system (2); taking a = l for simplicity, we 
present system (2) in the form: 



{x + y)-{x + yf-b{x + y)z + b{px + sy)z = 0, 
-d +bM{x + y) -b{px + sy){M+l) =Q, 
y + lx-y{x + y)-b{l-s + es)yz = 0, 
p = p{z) = {l-s)exp{-kz) + s 



(6) 



It follows from the first two equations of (6) that 

0. 



,2 , , {l + M-bz) dz 
(x + yf-{x + yy 



M+1 M+1 
The solution to the last equation is 

„ I ^^_ l+M-bz+^/[,l+bz-Mf-4d(l+M)z i+M-bz c„ 

y — 2(i+M) - 1+M ■ 

0<Xe-^y^< \ and Ze < — r — . 



(7) 
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It means that all possible non-trivial equilibria of the 
model are placed in a bounded area of the variable values; 
the amount of hosts is comparatively small (<1) and the 
amount of viruses is restricted by the virus reproduction 
rate M. 

Solving system (6) we find coordinates of non-trivial 
equilibrium B(Xe,ye, Ze) such that x = Xe,y = ye are expressed 
via z = Ze 

-d + bf^{z){M{l-s)-s) 



b{l+M){p-s) 
d-bf^{z){M{l-p)-p) 



, where 



b{l+M){p-s) 

1 + M-bz ±\J{'1+ M-bzf-4:ad{l + M)z 
2a{l+M) 



(8) 



and z -coordinate solves the equation: 

{-1 + l + af^) {d-bf^{M{l-s)-s)) 

+{des-b^f^{l-p){M{l-s)-s) + h{d{\-p) (9) 
-ef^{M{\-p)-p)s))z=Q. 

Analysis of the system (1) with a >0, (3) showed that 
there exists such threshold Ivf^ (depending on the model 
parameters) that the equilibrium B is stable if M < Nf'^; 
when M increases and intersects the threshold M = NF^, 
the equilibrium B loses stability and a stable limit cycle 
appears in the system. We summarize the results of this 
analysis in 

Statement 3. For a wide range of (fixed) parameters I, 
e,s,b,0<k< 1, model (1), (3) with a = 1 has a stable equi- 
librium for 0 < M < Ivf and stable oscillations for M > M", 
which appear due to supercritical Hopf bifurcation at 
M = M"'{1, e, s, b). 

See Methods for sketch of the proof 

Examples of the evolution of trajectories and phase 
portraits as the value of M increases are given in Figures 2, 
3, 4, 5. Figure 2 shows trajectories of the system that tend 
to a stable equilibrium when the parameter M is below 
the bifurcation threshold. Figure 3 demonstrates that the 
system arrives at a stable limit cycle when M intersects 
the bifurcation threshold, and Figure 4 shows that the 
established regime at the steady state is very close to 
periodic oscillations. 

Oscillations that are observed in Figure 4, with M 
above the bifurcation threshold IS/f', are very close to but 
not perfectly periodic. This peculiarity of the trajectories 
can be explained in the following way. The Hopf the- 
orem for a 3-dimension system states (see, e.g., [41], 
ch.5) that there exists such one-to-one transformation of 
the initial variables x — > x*, y — > /*, z — > z* that two of new 
variables (say, x'' and j*) show stable periodic oscillations 



but the third one does not and is governed by a separate 
equation. For this reason, the trajectories of initial 
variables, which are functions of x*, y", z*, may be non- 
periodic and may show not exactly periodic and even 
quasi- chaotic oscillations. 

It should be emphasized that, when the value of param- 
eter M crosses the bifurcation boundary, the qualitative 
behavior of the system changes sharply: as M increases, 
the behavior of the system becomes more and more 
"chaotic", with increasing non-regularity of the shapes 
of the trajectories (Figure 5). At large M, phase curves 
fill a surface in the [x,y,z) -space (see Figure 6, left 
panel, in contrast to Figure 4), and the phase portrait of 
the system reveals sharp, quasi-chaotic oscillations (see 
Figure 6, right panel). 

Discussion 

The bifurcation analysis of the virus-host system described 
here reveals complex, and in particular quasi-chaotic, 
oscillation regimes that so far have not been previously 
observed experimentally or in agent-based models of the 
CRISPR-mediated adaptive immunity [26,28-35]. The 
patchy distribution of the CRISPR-Cas systems across 
the microbial diversity, and even within relatively narrow 
groups of bacteria [10,42,43], implies complex dynamics 
of virus-host coevolution. Here we show that, in order to 
detect such non-trivial co-evolutionary regimes, the model 
has to be sufficiently complex, or put another way, less 
unrealistic than toy models that deal with a single host 
population and a simple dependence (or independence) of 
immunity on the amount of the virus, such as the two- 
component models examined here. The key factors for the 
appearance of quasi-chaotic oscillations are the non-linear 
dependence of the immunity on the amount of viruses 
and the three-dimensional phase space of the model 
which divides the hosts into the immune and suscep- 
tible populations, so that the system consists of three 
components (Table 1). In a conceptually similar manner, 
chaotic behavior has been observed in a model of 
bacteriophage-host interaction where the complexity of 
the system is introduced through inclusion of the time- 
delayed formation of consumable bacterial debris as a 
result of cell lysis [44,45]. 

It should be emphasized that the three-component 
models require the immune host to persist in the popula- 
tion, hence the results obtained here are only applicable to 
adaptive immunity that involves stable inheritance, i.e. the 
Lamarckian mode of evolution [13]. The CRISPR-Cas 
systems that function by introducing unique, directed 
modifications into the host genomes present the most 
straightforward case of such Lamarckian immunity [13,34]. 
Nevertheless, other adaptive immunity systems, in par- 
ticular the piRNA mechanism of transposon restriction 
in the animal germline, function on similar principles 
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[46]. Moreover, the siRNA branch of eukaryotic RNA 
interference, which is nearly universal among eukaryotes, 
also encompasses a substantial heritable component, 
albeit in this case via the epigenetic route [37,47-49]. 
The quasi-chaotic regime of virus-host coevolution de- 
scribed here might be relevant for these immunity systems 
as well. 

Both the Malthusian and the logistic versions of the 
three-component model demonstrate quasi-chaotic oscil- 
lations but they appear via distinct mechanisms (Table 1). 
The Malthusian version has no stable modes, with un- 
stable equilibria. For a wide area of the parameter values, 
the trajectories lie in a bounded domain, and for this rea- 
son, these trajectories demonstrate (quasi)chaotic behavior 
for all reasonable values of the virus birth rate M. In 
contrast, the logistic model has stable modes. There exists 
a critical value, namely a threshold of the virus birth rate 
Af" at fixed values of all other parameters, such that 
the system tends to a stable non-trivial equilibrium at 



M < M". When M > Nf\ the qualitative behavior of the 
system changes sharply. A limit cycle and (almost) 
periodic oscillations appear. With the further increase 
of M, the limit cycle turns into a "surface", and the behav- 
ior of the system becomes more and more "chaotic", with 
increasing non-regularity of the shapes of the trajectories. 
When M> > Af, the system transits to a regime of quasi- 
chaotic oscillations. 

Thus, the models described here make concrete, quan- 
titative predictions that can be directly tested using an 
experimental set-up for virus-host coevolution [27,50]. 

Conclusions 

Comparative genomics as well as direct experiments reveal 
notable evolutionary instability of the CRISPR-Cas systems 
[36,50-52]. It is common for closely related isolates of 
bacteria or archaea to differ with respect to the presence 
or absence of CRISPR-Cas, indicating that this immune 
system is easily lost and gained via HGT. Rearrangements 
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Figure 3 A limit cycle appears in the system (1), (3); M = 98.226 > M". 




of the CRISPR-Cas loci are also extremely widespread 
among microbes [43,53]. Given this evolutionary plasticity, 
complex and in particular quasi-chaotic regimes of virus- 
host coevolution revealed here appear to be plausible and 
potentially important for the evolutionary outcomes. In 
light of the current, rapidly increasing interest in CRISPR 
research, experimental validation of such regimes could 
be a realistic prospect. 

Methods 

Two-component model with p = const 

The host-virus dynamics may be considered within the 
framework of the Volterra-type models 

= x{l-ax-bz{l-p)) =P{x,z), 
z = z{-d-bxp + bMx{l-p)) = Q{x,z) 
x>0, z>0;a>0,b,d,M > 0 

(Ml) 

If the parameter p is constant and a = 0 then model 
(Ml) is Hamiltonian (conservative) with the Hamiltonian 
G{x, z) = In \z\ +d\n \x\ - b{M(l - p) - p)x - b{l - p)z. If 



Q < p < then the system has a saddle in the origin 
O and a center in the equilibrium 



b{M{l-p)-p) 



1-ax 
b{l-p) 



and B\x 



If a > 0 (logistic case), then for constant 0 < p < 
system (Ml) has a saddle in the origin, equilibria A{l/a, 0) 

d 2 _ b{M{l-p)-p)-ad \ . „ jljU i D 

b{M{i-p]-p) . 2 - b^(i.p)iM{i-p)-p) J ' equuiDrium a 
belongs to the positive quadrant {x, z) if b{M{l -p)-p)~ 
ad > 0. In this case, £ is a stable node/spiral and j4 is a sad- 
dle , A is a stable node if B is not positive (see, e.g., [54]). 

Two-component model with p = p(z) 

Malthusian version of the model (Ml) (fl = 0). 

This system has equilibrium in the origin (0,0), which 
is a saddle; z -coordinate of any other equilibrium (i,z) 
has to be a root of the equation 1 - bz{l - p{z)) = 0 and 
^ = biM-p{-z){i+M]) ■ The point {x,z) is positive only if 
p{z) <M/{1+M). 
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Proposition 1. Ifp^iz) < 0, then the Malthusian model 
has only one non-trivial equilibrium B(x,z), which is an 
unstable node/spiral for every parameter values. 

Proof. The Jacobian of the system in the equilibrium 



{X,Z) IS 



0 hx{-l + p{z) + zp^{z)) 

bz{-M + {I + M)p{z)) -b{l+M)xzp^{z) 



]{x,z) - 

and its determinant and trace are: 



Det(J(x,z)) = b^xziM-il + M)p(z))(l-p(z)-zp^(z)), 
Trace{J{x,z)) = -^^(1 + M)xzp^{z). 

(M2) 

If the function p{z} decreases monotonically, i.e. pj^z) < 0, 
then p(z) and monotonically increasing function h{z) = 1- 
^ can intersect only once. Thus, equation 1 - bz{l -p(z)) = 0 
has only one root z , and the Malthusian model has only 
one non-trivial equilibrium Next, Det{J{x,z)) > 0, 

Tr {J{x, z)) > 0 for positive x, z if Pz{z) < 0. Thus B{x:,z) is 
unstable node or unstable spiral. 

Logistic version of model (Ml) (a = 1). 



Equilibria of system (Ml) with a = l and p = p{z) defined 
by (3) are the points (0,0), A(1,0), and B{x,z) where coor- 
dinates X, z solve the system 



d 



b{M-p{z){l^ 
p{z) <M/{M+l) 



M)) 



l-x-bz{l-p{z)) = 0, (M3) 



Denote ^(z)= ^r^^Tsil) > *en z is a root of the equa- 
tion p{z) = h{z). Solutions of this equation are the points 
of intersection of the curves p{z) and h{z); up to two 
equilibrium points Bi{x\.,Z\).,B2{x2,Z2) can appear in 
the model with parameter variations. Denote J{x, z) the 
Jacobian of system (Ml), (3) with a = 1: 



J(x,z) = 

7(0,0) : 



l-2x-bz{l-p{z)) -bxil-p{z)-zp^{z)) 
-bz{-M + (1 + M)p{z)) -d + bMx(l-p(z))-bxp(z)-b(l + Mjxzpjz) 



1 0 

0 ~d 



-h{l-p{0)) 
-d + bM(l-p(0))-bp(0) 



0 

-d-i 



Thus 0(0,0) is a saddle and A{1,0) is a stable node for 
all parameter values. 
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z 




Figure 6 Left panel: phase curves fill a surface, M = 500; right panel: total host population against the virus population, 1 000 < t < 1 00000. 
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Table 1 The analyzed models of the coevolution between viruses and CRISPR-Cas-carrying hosts 







Two-component model 








Malthusian 


Logistic 


p 


= const 


Center M > p/O-p) Damp 


ed oscillations M > p/(1-p) 


p 


= p(z) 


No stable equilibrium z > 0 Damp 


ed oscillations M >A/I* (finite or infinite basin of attraction) 






Three-component model 








Malthusian 


Logistic 


p 


= const 


Damped oscillations M > p/(1-p) or s/(l-5) <M < p/(l-p) and / > /(e) Stable 


equilibrium for M > 


p 


= p(z) 


Quasi-chaotic oscillations at all M Damp 

quasi- 


ed oscillations at /M < M^''; periodic oscillations at M > A/I^'; 
chaotic oscillations at /M > > M^^ 



Now let us consider the system consisting of equation 
(M3) supplemented with the equation Det{J{x,z)) = 0. 
In the case Ai = Trace{J{x,z))*0 this system defines co- 
ordinates of saddle-node point B{x^z) whose second 
eigenvalue A2 = 0. The last system supplemented with 
the equation Tmce{J{x,z)) = 0 defines an additional 
degeneration in the system in the point B{x,z) where 
Ai=A2 = 0. 

Lemma 1. For a wide range of fixed parameters b, d, s 
there exist a point [M, /c) in the {M, k) -parameter plane 
such that the Bogdanov-Takens bifurcation of co-dimension 
2 is realized in the model (Ml), a = 1 under variations ofM 
and k close to (M,/c). The values (M^k) and coordinates 
of the equilibrium B{x^z) where x = x{M,k),z = z 
(M,/c) are defined by the system consisting of equations 
(M3) and equations Det{J{x.,z)) = 0, Trace{J{x,z)) = 0 
(see, e.g , [41];. 

We have found this bifurcation for some reasonable 
fixed values of the parameters d, b, s with the help of 
computer package LOCBIF [55]. Using the Lemma and 
Proposition 1 we prove the following statement. 

Proposition 2. (1) System (Ml), (3) has equilibria: the 
saddle 0(0,0) and stable-node A{\, 0) for all positive 
values of the parameters b, d, M, k = 1, s; 

(2) there exists positive M* such that 

a) for M < M* the system has only the equilibria O and 
stable equilibrium A; 

b) for M > M* the system has two more equilibria, a 
saddle Bi{xi,Zi) and a stable topological node/spiral 
£2(^2,22); 

c) there exist M > M* such that for M> M the spiral 
£2(^2,^2) is placed inside an unstable limit cycle. 

The bifurcation diagram of the system under variation 
of the parameter M is shown in Additional file 1. 

Three-component model: proof of Statement 1 

Let us formulate Statement 1 in more details: 



{1) Equilibrium 0{x = j = z = 0) has eigenvalues 
Ai(0) = 1, A2(0) = - d,X-i{0) = 1 - /; it is unstable 
for all parameter values; 

(2) If a > 0 then equilibrium A(0,l/a,0) has eigenvalues 
AM) = - I, HA) = - 1, X3(A) = -"'^+^-^(1-^)-^) ; it is 
stable for M < and unstable for M > . 

Proof. Jacobian of system (1) is of the form: J{x,y,z) = 
{aij)i,i = \, 2, 3, where 

= 1-1- ax- a(x +y) - bz{l - p{z)), ai_2 = - ax + esz, 
«i,3 - esy - bx{l - p{z)) + bxzpj^z); «2,i = 1- ay, 122,2 = 1-ay- 
a{x +y)-b{l- s)z- esz, <22,3 = -b{\- s)y - esy, 

^3,1 = bMz{l-p{z))-bzp{z)Tai,2 = bM{l-s)z-bsz, 

fl33 = -d + b(M{x + y)-(M + l){p{z)x + sy))-bxzpjz)(l + M). 

Pz{z) = 0, if = const and Pz{z) = l<{s - p{z)), if p{z) = (1 - s) 
e-''' + s. 

/!-/ 0 0 \ 

7(0,0,0)= / 1 0 , 
\ 0 0 -d J 

/ -/ 0 bes/a \ 

7(0, l/ffl,0)= -1 + / -1 b{-l + s-es)/a 

V 0 0 -{ad-b(M(l-s)-s)/a/ 

So, the equilibrium 0(0, 0, 0) is unstable for both 
Malthusian and logistic models, the equilibrium ^4(0, ^ , O) 
of logistic model is stable for M < and unstable for 
M > j^Y^ . The Statement is proven. 

Three-component Malthusian model with constant 
immunity p, p>s>Q 

Let us analyze the dynamics of model (1) depending on 
the parameters /, e. The system has trivial unstable 
equilibrium 0(0,0,0). Let firstly p = s. Then system (1) 
by changing of variables u = x + y, z = z is reduced to 2- 
component Malthusian system (Ml) with respect to u, 
z - variables. For s = p < the orbits {u, z} of this 
system belongs to the positive quadrant, the non-trivial 
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equilibrium u = , z = ^^is] ^ center, i.e., it is 

located inside closed orbits (similar to model (Ml)) and 
trajectories demonstrate periodic oscillations. In variables 
x{t),y(t), z{t) model (1) also demonstrates periodic oscUla- 
tions; the model has equilibrium B{Xe, je, Zg) where 



b{M{\-s)-s){bl{\-s) + esy 

dl{\-s) 
(M(l-s)-s)(W(l-5) + e5)' 
1 



(M4) 



h{\-s) 



If p> s then any non-trivial equilibrium (Xe, yet Ze) of 
model (1) is such that the coordinate solves the quad- 
ratic equation 

(1-/) + b{-2 + 1 + P + s-es{l-l))z + b^{l-p){l-s + es)z^ = 0, 

(M5) 

and coordinates Xi.,ye can be expressed via z = z^, as follows: 

-d{l-b{l-s)z) _ d{l-b{l-p)z) 
b{p-s){l + M-bz)' b{p-s){l+M-bz)' 



Let 1(e) = (-M + P + Mp)(i + M j.^tfli) ■ This equation defines 
a boundary line F = (e, /(e)) in the parametric plane (e, /). 

Proposition 3. 

1) System (1) with p = const and a = 0 has at most one 
positive equilibrium B{Xe,ye,z^ where 

defined by formulas (M5), (M6) ion s <p and by 
formula (M4) for s = p; 

2) the system has a single positive equilibrium B (Xe,ye, z^) 
if and only if one of the following conditions holds: 



^)s<p<^^; 
h)s<j^^<pandl>l{ey, 



3) in the case a) the equilibrium B is asymptotically stable; 

4) for s =p the system demonstrates periodic oscillations 
of variables x{t),y{t) while z(t) tends to a stable 
value for a wide domain of initial values {xqJo^Zo) 
close to the equilibrium B. 

Proof. 

Taking the solutions of quadratic equation (M5) in the 
form 



z = zi(/,e) = 
z = Z2(/, e) = 



2-l-(p + s) + es{l-l) + Vd 
2b{l-p){l-s + es) ' 

2-l-{p + s) + £5(1-/) + Vd 
2b{l-p)(l-s + es) 



where D = P(l - sf + {p - s + esf - 2l(p{l -s + 2es) - s{l + 
e - 5 + 65) we can easily verify that both "branches" Zi(/, e), 
Z2(/, e) are real for any positive (e, t) because the expression 
under the radical is non-negative. The branches Zi(/, e),Z2 
(/, e) are positive both if / < 1 and only Zi(/, e) is positive if 
/>1. 

Analysis of formulas (M5), (M6) shows that only the 
branch Zi(/, e) can define positive coordinates of the 
equilibrium x^ = x(zi),ye = yeizi). Substituting Zi(/, e) into 
formulas (M6) we obtain that Xe{e, l),ye(e, I) are positive 
if the point (e, /) in the parametric plane is placed above 
the boundary line F given by equation 

(-1 + p){l + (-1 + e)s){l{M{-l + s) + s) 

+ {M{-1 +p)+ p){{-l + e)s + M(l + (-1 + e)s)) = 0. 

Statements 1 and 2 of the Proposition are proven. 

Let us analyze a stability of equilibrium B {Xi,,ye>z^ of 
the system. For p = s characteristic polynomial of the 
system in the point B, whose coordinates are given 
by (M4), is of the form E{hq) = Det{{J{B)-piQl) = 



b(l-s) 



. Thus, two eigenvalues of the point 
(M6) 2j.g imaginary, /^g^ '^ = ±i\/~d. and the third is negative. 



1^0^ = ~ '"'X^[i-s) • Thus, statement 4 is proven. 

We show now that for p> s the point B is a sink, i.e., 
its eigenvalues have negative real parts (more exactly, one 
eigenvalue is real negative, and two others are complex 
with negative real part). 

Introduce the parameter a = p - s and write the right 
hands of (1), a = 0 in the form: 

P{x,y,z) = x-lx-b{l-a-s)xz + esyz, 

Q{x,y,z) = Ix + y-b{l-s)yz-esyz, 

R{x,y,z) = z{-d + bM{{l-a-s)x + {l-s)y) 
-b{{a + s)x + sy)). 

From the condition R{x,y, z) = 0, Q{x,y, z) = 0 we can 
express x = Xe,y = ye via z = z^: 

-{d{l-b{l-s)z + esz) 
^ b(a(\ Jr M){\-b{\-s)z ^ esz) + (M(l-s)-s)(-l + / + esz + fe(l-s)) ' 
dl 

^' ^ b(a{\ ^ M)(\-b{\-s)z ^ esz) + (*f(l-s) +s)(-l + / + esz + fe(l-s)) 

Substituting x = Xe,y = ye to P{x, y, z) we obtain: 



P{x„y„z) ^ H{z) = 



d{lll-b(l-s)z)-[l-Hl-s)z-esz}ll-b{l-s-ix)z)) 
b{iz{l + M){l-b{l-s)z + esz) + (M(l-s) + s)(-l + / + esz + te(l-s)). 



Solving the equation H(z) = 0 we get two roots 
Zei, Ze2', supposing z = Zo + ah^ and expanding H{z) in 
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series by a we get the two roots up to o{a) in the 
form: 

Zel = ZOI + ahiz, Ze2 = Zo2 + 0Ch2z, 

1 , es 



■i(M(l+ii)-s)+t.(l-s)(i(M(l+ii)-a)-rf(l+M)(l-a)) (l-i 



Zoi 



Z02 



b{l-s) 

l-l 



h 



ho. 



^(1-5)^(^(1-5)+ 65)' 
W(l-/) 



^(1-5) +es' {b{l-s) + es){bl{l-s) + es)' 

(M7) 

Notice, that the root z^i is positive, whereas 2^2 is positive 
only for 0 < / < 1. The coordinates x = xjfic), y = yXa) of B 
are both positive only for z = z^- 

Letting x^x = Xqi + ahi^i Jei = Joi + och^y, we get 



des 



Xoi - 



b(M(l-s)-s)(bl{l-s) + es)' 



^ _ -des{e^s^(l + M) + b^el{l-s)(M{l-s)-s) + bel{l + 2*f{l-s)s-2s^)) 
" " h{U{l-s)-sf(bl(l-s) + est ' 



dl(\-s) 
b(M(l-s)-s){blil-s) + es) )' 

_ dels{es-h{l-s){M-el{l+M)(l-s)+s + Ms)) 
" b(M(l-s)-sf{bl{l-s)+esf ■ 



(M8) 



Thus, coordinates of the equilibrium can be presented 
in the form {Xe,ye, Ze) = {x\e + och\x,y-Le + cihiy, Zie + a/ziz), see 
formulas (M7), (M8). 

We apply the method of small parameter expansion to 
verify stability of the equilibrium B^. Characteristic poly- 
nomial of the system can be presented in the form 



E{fi) ^ Det{{J{B,)-fiI) ■ 



bil-s) 



h{l-sy{M{l-s) + s)bl{l-s + es) ) 

where /j. = /4o + ah^^ and 

z{fi) = {-{b^l{d{l + h^{l-s)) + i4a{l^o{^-3h^{l-s)) 
+2lh^{l-s)){{{l-s)^M{l-s)-s) 
-es^{d + fio{fio + 2/!f,(l-s)))((l-5)2(M(l-s)-s) 
+be{l-s)s{fig{fif,{-l-3h^{l-s)) 
-m^{l-s)){M{-l + s) + s) 
+d({-l-h^{l-s)) (M(-l + s) + s) + 1{M{-1 + s) 
-^o(l + A^)(-l + s)+5))))). 

Let now fi^ = -d. Substituting this value to Z(/<) and 
solving equation Z(/^) = 0, we find 

^ _ (bdels(M + fif,(l + M)(-l + s)-s-Ms)) 

~ (2{Af(-l + s) + s)(&/(-l + s)-es){&(rf-//<o){-l + «) + es/^o)) 



< 0 



and Re(/,„)=5^( 
for reasonable parameter values. Thus, equilibrium B^ 
is asymptotically stable for s <p. 
The Proposition is proven. 

Three-component logistic model with constant immunity p 

The logistic version of model (1) with p = s is reduced to 
2-component logistic system (Ml) with respect to variables 
u = x + y,z. For s= p < it has two equilibria, A(0, 1/a) 

" " - 'T'lrMn'T'. ) • According to 
Proposition 2, these equilibria are stable in different 
parameter domains. 

In variables x,y, z the equilibrium £(Xe, je> ^e) of system 
(1) has coordinates 

des{b(M{\-s)-s)-ad) 
~ b{M{\-s)-s){b^l{\-s)(M(\-s)-s) + es{b{M{\-s)-s)-ad)) ' 
_ bdl{l-s) 

~ b{M{l-s)-s)(b'^l{l-s)(M(\-s)-s) + es{b{M{l-s)-s)-ad)) ' 
_ b{M{l-s)-s)-ad 



and £|M = ^(^(j-_^)_^) ,z 



r(l-s)(M(l-s)-s) ' 



(M9) 



Let phe a constant, p>s>0 . According to Statement 1, 
the system has trivial equilibrium 0(x = 0,^ = 0,2 = 0), 
which is unstable for all parameter values, and the equilib- 
rium A[x = 0,y = l,z = 0), which is stable if Al < 

and unstable if > fpifj ■ The system has also a non- 
trivial equilibrium B whose x,y - coordinates are expressed 
via z -coordinate: 

-d + bf^{z){M{l-s)-s) 



y = - 



b{l+M){p-s) 
d-bf^(z)(M(l-p)-p) 



b(\+M)(p-s) 



where 



^ 1 + M-bz ± \J (1 + M-bz) -^ad(l + M)z 
^ 2a(l + M) 

(MIO) 

and z -coordinate solves the equation: 

(-1 + / + af^) {d-bf^{M{l-s)-s)) + {des-b^f^{l-p){M{l-s)-s) 
+b{d{l-p)-ef^{M{l-p)-p)s))z = 0. 

(Mil) 

Denote 



(d(J-l+af^) +bf^(M{l-af"-l{l-s)) + Is)) + (d-bf^M)(b(l-s) + es)z 
bf^{l+M)(l-af^-(b(l-s) + es)z) 



then z is a root of the equations p = h*{z). Solutions of 
the latter equation are the points of intersection of the 
line 0 <p <1 and the curves /? (z). Two cases with 
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different values of the parameter M are presented in 
Additional file 2. It demonstrates that at most two positive 
values of z are possible and hence the system may have at 
most two positive equilibria, whose x, y - coordinates are 
expressed via z by the equations (M9). 
Let us define the critical value of the parameter M, 

j^tc ^ ad±bs 
A(l-s) 

Proposition 4. For s<p < and fixed positive values 
of d,b < 1, /, e system (1) with a =1 has a single stable 
equilibrium A{x = Q,y = ^,z = 0) if M < Ivf'^ and a sin- 
gle stable equilibrium B[xQ,yQ, Zq) if M > M^'^; here the co- 
ordinates {xci,yo, Zq) are defined by formulas (MIO, Mil) 
if s < p and by (M9) if s = p. Transcritical bifurcation 
"changing of stability" of the points A and B happens as 
M = 

For p = s system (1) has nontrivial point B, whose 
coordinates are expressed by (M9); it is easily to 



b{M{l-s)-s)-ad 
b'^(l-s)(M(l-s)-s) ' 



verify that Ue = Xe+y, = i,^M(ts)-s) ' ■ 
and 1 - a{Xe + - b{l - s)Ze = 0. Hence, z^, > 0 if > M^"; at 
this condition the equilibrium A looses stability according 
to Statement 1. 

Characteristic equation at the equilibrium B for p = s is 
of the form: 

Det{{J{B)-fiI) = -(-1 + l + fi + a(x, +yj + b{l-s)Ze + sz,)* 
(f? + b'^(l-s)(M(l-s)-s)(x, + y^)ze + + 2a(xi,+y^) + fee(l-s))) = 0. 

Accounting the equality 1 - a{Xe+ye) - b{\ - s)Ze = 0, we 
obtain from the first term of the last equation: = 1 - / - 
aUe - b(l - s)Ze - esZe = - / - 5Ze < 0, and from the second 



term pi^ 3 



-au ± •fj -i{b{M{l-s)-s)-ad)Ue + (au)' 



For 



M'^ f4.2 = 0, ^3 = - AM < 0. Hence, the transcritical (pitch- 
fork) bifurcation happens with points A, B [41], ch.7. 
When M > M'^ the eigenvalues ^2,3 have negative real 
parts. So, the equilibrium point B is stable for p = s and 
> M*"^. Due to continuity arguments the point B is 
also stable for /7 > 5 if is close to s. A general case p>s 
was verified by computation of eigenvalues of complete 
characteristic equation in B. 



p{z) — > 5 < 1 fls z — >■ 00, so that two general cases of mutual 
placing of p{z) and h{z) for / > 1 and / < 1 are possible, 
which are shown in Additional file 3. Equation (M12) 
has up to two ("small") positive roots z = Ze if / < 1 (see 
Additional file 3a) and one ("large") positive root z^ if / > 1 
(see Additional file 3b). Equilibrium values Xe{z),ye(z) 
defined by the formulas (M6) have different signs for 
small values of Ze and are both positive for large Ze- Ac- 
cordingly, the system can have only one positive equilib- 
rium. Eigenvalues of this equilibrium can be computed 
from the equation Det{J{Xe,ye,Ze) - fil) = 0. Using the 
LOCBIF software [55], we show that one eigenvalue is 
real and negative but two other are complex with a posi- 
tive real part. Thus, this equilibrium is unstable. Statement 
2 is proven. 

Three-component logistic model with the immunity p = p(z) 

Proof of Statement 3 

Let BM{x,y,z) be a non-trivial stable equilibrium of 
model (1), (3) whose coordinates {x(M), y{M), z(M)) 
depend on the parameter M and satisfy system (2); let 
P{pi) = Det{J(B) - ^7) = 0 be the characteristic polynomial of 
the system around Bm- The supercritical Hopf bifurcation, 
corresponding to changing of stability of Bm accompanied 
by appearance a stable limit cycle happens in the system 
when for some M a pair of eigenvalues becomes imaginary 
and certain conditions of non-degeneracy are fulfilled (see, 
for example, [41]). The Hopf bifurcation is supercritical if 
the first Lyapunov value becomes negative. The existence 
of this bifurcation in logistic version of model (1), (3) was 
verified using LOCBIF [55]. The program numerically finds 
coordinates of equilibrium with imaginary eigenvalues 
under variation of parameter M and one more parameter 
of the model (e.g., e, /, or s) for fixed values of other pa- 
rameters, checks the sign of the first Lyapunov value and 
verifies the conditions of non-degeneracy formulated in 
the Hopf theorem. Applying this software we have found 
the parameter curves of Hopf bifurcation e//(M), Ij-^M); Sh 
{M) for fixed values of other parameters of the model 
(see Additional file 4); it proves the assertions of the 
Statement. 



Three-component Malthusian model with the immunity 
p = p(z); proof of Statement 2 

The coordinates of non-trivial equilibria are given by for- 
mulas (M5) and (M6). Solving equation (M5) with respect 
to p = p{z), we obtain the equation for z-coordinate: 

'^^ ' ^ ' te(l-(fe(l-s) +es)z) 



=-Hz) 



(M12) 



Evidently, the z-coordinate of possible equilibrium 
does not depend on the parameter M. Next, h{z) 1, 



Reviewers' reports 

Reviewer 1: Sandor Pongor 

Comment: The CRISPR-Cas system is of high theoretical 
and practical interest. On the practical side it is important 
for designing genome engineering tools, on the theoretical 
side, it provides a noteworthy example of Lamarckian evo- 
lution. By inserting virus-derived spacers into CRISPR re- 
peat cassettes, microbes preserve in their genomes the 
signature of viruses that attack them, which in turn they 
pass on to their offspring. The evolution of this system is 
practically intriguing and has been the subject of several 
agent based modeling studies. As the authors point out. 
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agent systems can be used to model various levels of 
complexity, however they do not permit a full and rigor- 
ous mathematical analysis of all possible behaviors. 
Berezovskaya and associates present here a differential 
equation based model of the evolution of CRISPR-Cas 
based immunity. They use Lotka-Volterra type models 
that allow the analysis of Malthusian and logistic regimes 
and show that the models display complex quasi-chaotic 
oscillations. Such complex behaviors have not been found 
either by experiment or by the previous agent based 
mutations. The results are well underpinned and clearly 
discussed. The results are especially interesting example 
of how unexpected complexity emerges in a seemingly 
simple system. The text is very well written however the 
authors may want to check it for minor typos, e.g.: 

Page 3 challenge (a virus), "in contrast as opposed" to the 
random, undirected mutations in the Darwinian evolution- 
ary framework - one of them is superfluous. 

Page 4 after eqn 1 host reproduction is Malthusian if 
a = 0 or logistic if a > 0, and in both cases. - sounds 
like an unfinished sentence. 

Authors' response: We appreciate these comments. The 
proposed corrections have been made. 

Reviewer 2: Sergei Maslov 

Comment: The manuscript addresses a very interesting 
topic of the emergence of chaotic oscillations in popula- 
tion dynamics of phages and their bacterial hosts. To my 
knowledge, this topic is relatively unexplored in the litera- 
ture except for "Bifurcation analysis of bacteria and bac- 
teriophage coexistence in the presence of bacterial debris" 
by Ira Aviram, Avinoam Rabinovitch [44] . Authors should 
cite this paper and compare and contrast their results to 
findings of Aviram et al. 

Authors' response: As far as we can see, the work of 
Aviram and Rabinovitch is relevant only in a very generic 
sense. In the revision, we cite this paper along with a more 
recent publication of the same authors, and comment on 
the emergence of complex behavior in these models. 

Comment: I found the paper very hard to read and 
understand. In its present form it is unnecessarily heavy 
on mathematical details and terminology and light on 
biological insights. I would strongly recommend a *near 
complete rewriting* of the text of the manuscript that 
would delegate unnecessary mathematical details and 
theorems to supplementary materials and explaining the 
biological consequences of main findings. For example, 
on page 5 authors refer to their model as "conserva- 
tive". It is not at all obvious to the majority of even so- 
phisticated computational biology readers that in a 
conservative system small deviations from the steady 
state solution do not decay back to the steady state 
but persist indefinitely. Whenever possible authors 
should avoid using mathematical jargon and explain in 



plain English what their parameters/assumptions mean 
biologically. 

Authors' response: As per this comment and analogous 
comments of Dr. Kimmel (see below), the entire description 
of two-dimensional (or two-component, using the modified 
terminology), which included most of the mathematical 
detail, was moved to the Methods. There, we considered it 
appropriate to formally mathematically define "conserva- 
tive systems" and other concepts that might be unfamiliar 
to non-specialist readers. In the main text, biological inter- 
pretations of the parameters and assumptions were pro- 
vided on several occasions. We believe that these changes 
made the paper considerably more straightforward, and 
we appreciate these suggestions of the reviewers. In general, 
however, one has to face the fact that this is a mathemat- 
ical biology paper. Further, we beg to disagree that this 
paper is "light on biological insights". We think that the re- 
vealed complex behavior of the co-evolving virus-host sys- 
tems is a potentially important biological instant, and the 
reviewers do not seem to disagree. What is somewhat lack- 
ing are direct connections to experimental results. Such 
relevant results could come, first, from quantitative mea- 
surements of virus diversity and CRISPR-Cas prevalence in 
various habitats, and second, from laboratory co-evolution 
experiments. We expect that such results are indeed forth- 
coming and hope that the present facilitates their inter- 
pretation but at this time, there is little data for direct 
comparison. 

Comment: Another example of the same point: on 
page 5 authors introduce two steady state solutions A 
and B but do not explain what they mean biologically: 
phages co-exist with bacteria in A but die off in B. And 
this is just one example of the lack of biological inter- 
pretation of mathematical results happening throughout 
the manuscript. 

Authors' response: This specific statement has been 
moved to the Methods as per the suggestion of Dr. Maslov 
and Dr. Kimmel. The interpretation of these steady solu- 
tions given by Dr. Maslov is quite correct and thus, we 
presume, is obvious from the presentation. Especially in 
the Methods section, further clarifications seem unneces- 
sary. On several other occasions (see below), however, we 
did include additional biological interpretations. 

Comment: On the same page authors cite earlier studies 
with empirical results confirming that the fraction of im- 
mune bacteria p directly depends on the phage population 
z and not on the bacterial population x. A more detailed 
discussion of what the empirical data actually say would 
be beneficial here. 

Author's response: This discussion has been expanded 
and made more specific. 

Comment: In particular, I don't expect p(T) to instantan- 
eously trace rapidly growing phage population z(T), which 
seems to be a prerequisite for chaotic behavior reported in 
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this manuscript. What would happen when a delay or time 
averaging is added to the model? That is to say, what would 
happen if p(T) is a function of z(T-t_delay) or (in a separate 
ariant of the model) p(T) is determined by the time- 
averaged value of z(T) over T:T-t_average time interval? I 
would be particularly interested to see if chaotic dynamics 
would disappear for large enough values of t_delay or 
t_average. What is a realistic value of t_average compared 
to a single generation of lytic phage growth? 

Authors' response: These are interesting and important 
questions that, however, are beyond the scope of the present 
paper. 

Comment: I also request a small yet important change 
in notation: authors repeatedly refer to three-dimensional 
version of their model which I understood as a model with 
3 spatial coordinates. Indeed, special inhomogeneity is 
known to play an important role in phage-bacterial inter- 
actions (see e.g. [56,57]. However, what authors meant is 
simply a three-species or three-component model with 
phages, immune hosts, and susceptible hosts. To avoid 
confusion the term "three-dimensional model" while 
mathematically correct needs to be changed to "three 
component model" throughout the manuscript. 

Authors' response: We adopted this change in the 
revision. 

Comment: Are figures of plots in Figure 1, 2, 3, 4, 5 
necessary? I personally did not learn anything from them. 
Figure 6, Additional files 1, 2, 3 and 4 nicely illustrate cha- 
otic dynamic of the system. Perhaps they can be collected 
as multiple panels of just one figure? 

Authors' response: We agree, figures 1, 2, 3, 4, 5 have 
been moved to the additional files. 

Additional comments made in the second round of review 

Comment: On the other hand, steady state and simple 
time-dependent solutions to Lotka-Volterra equations 
for bacterial-phage systems have been studied for quite 
some time. Some classic references such as [22,23] have 
been overlooked and need to be cited in the manuscript. 

Authors' response: We agree, these are relevant refer- 
ences that are cited in the revised version of the present 
article. 

Comment: Population cycles and fixed points in modi- 
fied Lotka-Volterra equations have been also considered 
in [58,59]. Would authors predictions of chaotic (or quasi- 
chaotic) behavior persist in these systems? 

Authors' response: We do not see how to directly link 
these models with our approach (at least not without 
additional, extensive analysis) and therefore currently 
cannot make such predictions. 

Comment: Spatial (see e.g. [56]) and temporal [57] 
inhomogeneity of the environment is known to play an 
important role in phage-bacterial interactions. How much 
would it affect chaotic dynamics. These are not idle 



questions since they go to the heart of the question of 
how generic is the chaotic behavior reported in the 
manuscript. If authors believe they are beyond the 
scope of the current paper, perhaps, these questions 
should be mentioned in the discussion section as model 
generalizations/modifications that need to be performed 
in future studies. 

Authors' response: we fully agree that these are among 
desirable generalizations of the present approach. It is 
another matter whether, with the addition of these non- 
homogeneities, the model remains tractable. On very 
general grounds, given that here we have shown that 
the pseudo-chaotic oscillations only emerge in a system 
with certain minimal complexity (distinguishing suscep- 
tible and immune hosts is essential), we would expect that 
such oscillations only become more prominent in even 
more complex models. However, this is obviously only a 
conjecture at this point, we cannot be confident before the 
actual analysis is done. 

Comment: Throughout the manuscript authors consider 
only virulent (lytic) phages. Would there be any interesting 
modification of predicted dynamical patterns for temperate 
phages? 

Authors' response: As such, temperate phages, by defin- 
ition, do not kill the host, and therefore, even if lysogeniza- 
tion is prevented by CRISPR-Cas, as indeed has been 
reported [60], this seems to be irrelevant for the modeling 
approach described here. The situation certainly changes 
when it comes to prophage induction, against which 
CRISPR-Cas protects as well [60]. This case does not 
appear to be distinguishable from lytic infection within the 
approximations of the model. 

Comment: I appreciated authors following my request 
and renaming two/ three dimensional model into two/three 
component model. However, on page 11 and several other 
places in the manuscript old notation is still being used. I 
recommend authors do global search for "2D", "3D", 
and " dimensional" in the manuscript. 

Authors' response: This has been taken care of. 

Reviewer 3: Marek Kimmel 

This paper addresses the issue of co-evolution of a virus 
an and the immune system of a host, taking into account 
the dynamics of the virus and two types of immune 
systems: susceptible and resistant. The model is in- 
spired by a type of immune response (CRISP-Cas) in 
archaea and some bacteria. The dynamics is summa- 
rized by a system of 3 nonlinear ordinary differential 
equations (ODEs). The system seems to exhibit various 
dynamical regimes including some that are chaotic. 
This, according to the authors, provides some analogy 
to the known examples of the CRISP-Cas system 
behavior. 

The paper should be reorganized before it is publishable. 
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Major issues 

Comment: 1. The paper is written in a way which malces 
understanding it very difficult. A large portion of the 
paper is devoted to models which are inadequate in 
that they do not include sensitive and resistant im- 
mune systems, are therefore limited to two ODEs and 
cannot exhibit complicated dynamics. To make the paper 
readable, it should proceed directly to the point. The aux- 
iliary models can be moved to an appendix. 

Authors' response: This reorganization of the manuscript 
has been implemented as suggested. 

Comment: 2. Dynamics of the really interesting 3-ODE 
models is explored mostly numerically, if I understand 
correctly. In my opinion, more illustrative material might 
be provided, using the space available after removal of the 
2-ODE models. 

Authors' response: We carefully considered this sugges- 
tion but found that the comparison of Figures 2, 3, 4, 5 
and 6 was highly illustrative of the results for the 3-ODE 
models. The transitions between the outcomes depending 
on the parameters was made fully explicit in the revised 
description. 

Comment: 3. I am missing a detailed discussion of 
how the model is similar to experimental observations. 
What type of data are available? Time series are prob- 
ably most desirable. Is there any information in the data 
regarding sensitivity to parameter change and initial 
conditions, which are characteristic of chaos? 

Authors' response: As pointed out in our response to 
Dr. Maslov's comments above, the experimental data are 
simply not ready for detailed comparison. We would be 
more than happy to analyze time series or cite an appro- 
priate analysis but such experiments belong in the future. 

Comment: 4. Finally, why would this quite generic 
mathematical description be specific to the CRISP-Cas 
systems? 

Authors' response: We never claimed that this de- 
scription was specific to CRISPR-Cas. It is inspired by 
CRISPR-Cas but is applicable to any system of adaptive 
immunity with sufficiently long term memory, as pointed 
out both in the abstract and in the Conclusions. 

Detailed remarks 

Comment: 1. Page 9. It seems x and y are each composed 
of a number of different "immuno-types" of host individ- 
uals. The structure will be continuous and not two-point as 
it is now. What will be the consequences for the dynamics? 
Will it not be more regular because of smoothing effect of 
continuity? 

Authors' response: To the best of our understanding x 
is just one type. However, y indeed can be represented as 
numerous "immuno-types" if immunity to different viruses 
is considered separately. This situation has been explored 
within the framework of agent-based models [28,32]. Under 
the analytic approach used here, continuous distribution 



of "immune-types" would inevitably make the model 
intractable. 

Comment: 2. Page 15. Computational analysis usually 
cannot "show" that an equilibrium is stable. It may be at 
best consistent with stability. 

Authors' response: The revised version of the manuscript 
includes a comprehensive test for stability using LOCBIF. 
Thus, "show" (which does not mean "proven") seems 
appropriate. 

3. "Proposition 4" does not seem to be mathematically 
demonstrated. So, it is a Conjecture. 

Authors response: A sketch of the proof is given in the 
revision. 
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